Attach packages:

library(tidyverse)
library(here)
library(janitor)
library(raster)
library(sf)
library(tmap)
library(tmaptools)
library(gstat) # variogram

Grand Canyon DEM

Read in some raster data (GeoTIFF DEM)

Grand Canyon GeoTIFF from USGS: https://pubs.usgs.gov/ds/121/grand/grand.html

# Read in the data with raster::raster()
gc_dem <- raster(here("data","gc_dem.tif"))

Go exploring (basic plot, check CRS, projection, bounds)

# Look at it with base plot()
plot(gc_dem)

# Check CRS & bounds:
gc_dem@crs # Shows CRS: WGS84
## CRS arguments:
##  +proj=utm +zone=12 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0
gc_dem@extent # Shows extent (bounds)...notice that these seem odd (not units)
## class      : Extent 
## xmin       : 229520 
## xmax       : 412580 
## ymin       : 3982160 
## ymax       : 4100630
wgs84 = "+proj=longlat +datum=WGS84 +ellps=WGS84 +no_defs" # Just have this ready to copy/paste

# Reproject
gc_reproj = projectRaster(gc_dem, crs = wgs84, method = "bilinear")

# Then check: aha, now degrees we're used to
gc_reproj@extent
## class      : Extent 
## xmin       : -114.0417 
## xmax       : -111.9681 
## ymin       : 35.94488 
## ymax       : 37.04945

Crop the raster to a smaller area:

bounds <- as(extent(-112.4, -112.0, 36.1, 36.3), 'SpatialPolygons')

# Make the projection for "bounds" the same as for "gc_reproj":
crs(bounds) <- crs(gc_reproj)

# Then crop gc_reproj by the new bounds polygon: 
gc_crop <- crop(gc_reproj, bounds)

# Look at it:
plot(gc_crop)

Want to resample? Use raster::aggregate() to create lower res (larger cell) rasters.

See ?aggregrate (default is mean, fact is number of cells in each direction; can set 2 if wanting x/y to differ for aggregation)

# Aggregate: 
gc_agg <- aggregate(gc_crop, fact = 10)

# Then look at it: 
plot(gc_agg)

Get safe and sound back in ggplot:

# First, convert to a data frame:
gc_df <- as.data.frame(gc_agg, xy = TRUE)

# That `xy = TRUE` is important: retains the lat/lon information! 

# View(gc_df)

ggplot(data = gc_df, aes(x = x, y = y)) +
  geom_raster(aes(fill = gc_dem)) +
  coord_quickmap() +
  theme_minimal() +
  scale_fill_gradientn(colors = c("purple",
                                  "magenta",
                                  "orange",
                                  "yellow", 
                                  "white")
                       )

Let’s say we know that in this region, a certain species will only grow between 1000 and 1500 ft elevation. Create a subset of gc_crop that includes that habitat:

# First, make a copy
gc_hab <- gc_crop

# Set any cells outside of 1000 - 1500 to NA
gc_hab[gc_hab > 1500 | gc_hab < 1000] <- NA

# Plot.
plot(gc_hab)

# Cool!

Let’s make a bit nicer map with tmap:

tmap_mode("view") # Set to interactive viewing

# Make tmap: 
tm_shape(gc_hab) +
  tm_raster(legend.show = FALSE, palette = "plasma")
# In console, run tmaptools::palette_explorer() to view a Shiny app with other palettes!

Kriging Kansas Rain

Read in the KS counties shapefile with read_sf:

ks_counties <- read_sf(here("data",
                            "ks_counties",
                            "ks_counties_shapefile.shp"))

# View(ks_counties)

Do some exploring:

# Base plot
plot(ks_counties)

# Check CRS:
st_crs(ks_counties) # hmmmm none...guess we should set one!
## Coordinate Reference System: NA
# Set to EPSG 4326 (WGS84 datum):
st_crs(ks_counties) <- 4326

# Now check again:
st_crs(ks_counties)
## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"
# And replot:
plot(ks_counties)

# That looks more like Kansas.

Now, read in the KS rainfall data (this was from a single rainy day in 2017):

ks_rain <- read_csv(here("data","ks_rain.csv")) %>% 
  clean_names()

But currently, R has no idea that these are spatial points. We’ll convert it using sf::st_as_sf():

ks_sf  <-  st_as_sf(ks_rain, coords = c("lon", "lat"), 
                 crs = 4326)

# View(ks_sf) (rainfall = amt)
plot(ks_sf)

# Or in ggplot:
ggplot() +
  geom_sf(data = ks_counties) +
  geom_sf(data = ks_sf, aes(color = amt,
              size = amt),
          show.legend = FALSE) +
  coord_sf() +
  scale_color_gradient(low = "yellow", high = "red") +
  theme_void()

Now let’s try to predict rainfall all over Kansas based on those rainfall values.

For kriging, we’re gonna switch over to sp functions. Which means we need to get our data to ‘Spatial’ format, instead of sf. Eventually, I think these will all be merged…

Make the rainfall data a Spatial Points data frame (so we can use gstat/sp functions:

ks_sp  <- as_Spatial(ks_sf)
# class(ks_sp)

Make a grid that we’ll krige over:

# bbox(ks_sp) to check bounding box of the spatial points
lat <- seq(37, 40, length.out = 200)
long <- seq(-94.6,-102, length.out = 200)

# Then make it into a grid: 
grid <- expand.grid(lon = long, lat = lat)
grid_sf <- st_as_sf(grid, coords = c("lon","lat"), crs = 4326)
grid_sp <- as_Spatial(grid_sf)

# Check out your amazing kriging-ready grid:
plot(grid_sp)

Make a variogram from the rainfall data:

# Create the variogram:
ks_vgm <- variogram(amt ~ 1, data = ks_sp)

# Look at it: 
plot(ks_vgm)

# Fit the variogram model using reasonable estimates for nugget, sill and range:
ks_vgm_fit <- fit.variogram(ks_vgm, model = vgm(nugget = 0.2, psill = 0.8, model = "Sph", range = 200))

# Plot them both together
plot(ks_vgm, ks_vgm_fit) # Cool! So what are the values

# Just FYI: there are other models (Gaussian, Exponential) - how do those line up? 
ks_vgm_gau <- fit.variogram(ks_vgm, model = vgm(nugget = 0.2, psill = 0.8, model = "Gau", range = 200))

plot(ks_vgm, ks_vgm_gau)

# You can check the sum of squares of residuals for each: 
attr(ks_vgm_fit, 'SSErr') # 0.00214 (and could compare to other models...)
## [1] 0.002138374
# We'll stick with the Spherical model: 
ks_vgm_fit # Nugget = 0.102, sill = 0.954, range = 235
##   model     psill    range
## 1   Nug 0.1021677   0.0000
## 2   Sph 0.9537363 235.1416

Now, kriging!

ks_krige <- krige(amt ~ 1, ks_sp, grid_sp, model=ks_vgm_fit)
## [using ordinary kriging]
spplot(ks_krige, "var1.pred")

Let’s get it back into a format we’re used to (data frames & sf objects):

# Make a data frame from kriged predictions:
ks_df <- data.frame(ks_krige@data["var1.pred"], ks_krige@data["var1.var"],
                    ks_krige@coords) %>% 
  rename(longitude = coords.x1,
         latitude = coords.x2)

# Convert to sf object: 
rain_sf <- st_as_sf(ks_df, coords = c("longitude","latitude"))
st_crs(rain_sf) <- 4326
# Get Kansas outline to crop: 
ks <- read_sf(dsn = here("data","states"), 
              layer = "cb_2017_us_state_20m") %>% 
  dplyr::select(NAME) %>% 
  filter(NAME == "Kansas") %>% 
  st_transform(crs = 4326)

plot(ks)

# Find the intersection of the two: 
rain_sf_ks <- st_intersection(rain_sf, ks)

ggplot(rain_sf_ks) +
  geom_sf(aes(color = var1.pred)) +
  scale_color_gradientn(colors = c("white","yellow","magenta","purple")) +
  theme_minimal()